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ABSTRACT 

This  report  describes  a computer  program,  SLAM,  which  uses  ship 
slamming-impact  forces  and  structural  parameters  of  the  hull  to  calculate 
the  vertical  hull-girder  vibratory  response  in  terms  of  displacements,  accel- 
erations, bending  moments,  shear  forces,  and  stresses.  The  normal  mode 
method  is  used  so  that  the  user  can  calculate  only  the  modes  of  interest 
and  can  eliminate  the  rigid  body  and  higher  modes.  Modal  responses  are 
calculated  with  a time-marching  technique.  The  program  uses  a finite- 
element  model  of  a beam  which  is  suitable  for  conventional  monohulls. 

Modifications  to  accommodate  other  model  configurations  can  be  made. 

Test  problems  show  good  agreement  with  exact  solutions  for  a 
uniform  Euler  beam.  Sample  calculations  are  made  on  a MARINER- 
Class  hull  using  3 modes  and  again  using  10  modes. 

ADMINISTRATIVE  INFORMATION 

The  work  described  in  this  report  was  authorized  and  funded  by  the  Naval  Sea  Systems 
Command  (SEA  03412)  under  Project  S-F354  21,  Task  Area  S-F354  21  007,  Program  Element 
62512N.  Work  Units  1-1506-006  and  1-1568-102. 

INTRODUCTION 


BACKGROUND 

The  Naval  Ship  Research  and  Development  Center  (the  Center)  has  developed  a com- 
puter program  for  evaluating  the  impact  loads  and  the  main  hull  girder  response  associated 
with  slamming  of  a ship  at  sea  so  that  a more  rational  design  of  high-speed  naval  ships  can  be 
achieved.  The  program  will  calculate  impact  forces  as  a function  of  the  configuration  of  the 
ship  bottom,  speed,  draft,  trim,  and  state  of  sea.  Impact  forces  and  structural  parameters 
will  then  be  used  to  obtain  the  hull  girder  response  in  terms  of  relative  displacements,  accel- 
erations, bending  moments,  shear  forces,  and  stresses.  At  present  the  computer  program  is  in 
two  parts;  one  calculates  the  forces  and  the  other  the  response.  They  could  easily  be  com- 
bined. Results  from  the  overall  study  are  presented  in  Reference  1. 

This  report  presents  details  of  the  program  dealing  with  the  hull  response.  Hull  struc- 
tural characteristics  and  impact  forces  are  the  input  to  the  program. 

A review  of  experimental  data  about  slamming  shows  that  the  hull  response  is  most 
significant  in  the  lower  modes  of  vibration.  Also,  it  has  been  observed  that  although  the 


'ochi.  M.  K.  and  L.  E.  Motter,  “Prediction  of  Slamming:  Characteristics  and  Hull  Responses  for  Ship  Design,"  Society  of 
Naval  Architects  and  Marine  Engineers  Transactions  (1973).  A complete  listing  of  references  is  given  on  page  41. 
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forces  are  dependent  on  rigid  body  motion,  the  rigid  body  motion  itself  is  affected  very  little 
by  slamming. 

To  determine  whether  the  inclusion  of  shear  rigidity  and  rotary  inertia  would  increase 
the  accuracy  of  the  program  significantly,  Center  reports  about  calculated  natural  frequencies- 
of  ship  hulls  were  reviewed.  Three  were  found  that  contained  useful  comparisons.2-4  A sum- 
mary of  these  results  and  some  unpublished  data  are  given  in  Table  1 . 

Rotary  inertia  is  not  now  normally  included  in  vibration  calculations  at  the  Center,  and 
the  example  of  the  T-AGM-19  indicates  that  the  practice  is  justified.  Shear  rigidity  is  nor- 
mally included  because  the  calculations  usually  involve  at  least  the  first  five  modes  of  vibration, 
and  in  the  higher  modes  shear  rigidity  significantly  influences  the  modal  frequencies.  Another 
trend  which  is  reflected  in  the  table,  and  which  is  to  be  expected,  is  that  the  shear  rigidity  has 
less  influence  on  long  slender  ships  than  on  short  deep  ships.  Slamming  is  more  of  a problem 
on  long  slender  hulls.  To  make  the  program  as  general  as  possible,  however,  shear  rigidity 
was  included. 

It  is  recognized  that  the  accuracy  of  natural  frequency  calculations  is  only  one  factor  in 
the  accuracy  of  hull  response  calculations,  but  it  is  one  of  the  few  indicators  available. 

Many  of  the  requirements  of  this  program  are  similar  to  those  for  the  author’s  thesis5 
at  Catholic  University,  and  many  of  the  same  techniques  are  used. 

APPROACH 

Calculation  of  steady-state  hull  vibration  at  the  Center  has  normally  involved  lumping 
the  hull  parameters.  However,  the  advantages  of  a finite-element  model  are  enough  to  war- 
rant a different  breakdown  of  parameters.  A finite-element  beam  model  is  used. 

The  normal  mode  method  seems  particularly  suited  to  the  problem  of  ship  slamming  and 
is  used  in  the  program.  It  enables  the  user  to  calculate  only  the  modes  that  deserve  consid- 
eration and  to  eliminate  the  higher  modes  and  rigid  body  motion.  An  alternative  method 
would  involve  representing  the  buoyancy  of  the  ship  so  as  to  restrain  rigid  body  pitch  and 
heave. 

To  use  the  normal  mode  approach  it  is  necessary  to  calculate  all  of  the  natural  fre- 
quencies and  mode  shapes  of  interest.  The  next  major  step  is  to  transform  the  coordinates 


'Robinson,  D.  C\,  "Calculated  Natural  f requencies  and  Normal  Modes  of  the  Guided  Missile  Ouiset  USS  LONG  BEACH 
(C'G(N>-9)."  David  Taylor  Model  Basin  Report  2100  (Jan  1966V 

^Perkins.  R.  L.,  "Calculated  Natural  f requencies  and  Notmal  Modes  of  Vibration  on  Range  Instrumentation  Ship 
(T-AGM-19)."  David  Taylor  Model  Basin  Report  1997  (Jun  1965). 

4McGoldiick,  R.  T.  and  V.  L.  Russo.  “Hull  Vibration  Investigation  on  SS  GOPHER  MARINER,"  David  Taylor  Model 
Basin  Report  1060  (Jul  1956). 

5Antontdes,  G.  P„  "A  Computer  Program  for  Normal  Mode  Solutions  in  Structural  Dynamics,"  Master's  Thesis  at 
Catholic  University.  Washington.  D.  C.  (Dec  1970). 


to  obtain  the  uncoupled  equations  of  motion.  Each  equation  then  represents  a mode  of 
vibration.  The  modal  equations  are  solved  for  the  modal  displacements  in  the  modes  of 
interest,  and  these  are  transformed  back  into  physical  displacements.  The  accelerations,  bend- 
ing moments,  shear  forces,  and  stresses  can  then  be  found  from  the  displacements  and  ship 
parameters. 

For  the  solution  of  modal  equations,  the  program  uses  a time-marching  method  involv- 
ing finite  differences  with  respect  to  time.  The  modal  displacements  at  any  instant  are  calcu- 
lated as  a function  of  the  displacements  at  earlier  times,  system  parameters,  and  applied 
forces.  To  minimize  the  use  of  computer  storage,  all  of  the  required  quantities  are  calculated 
for  one  instant  of  time  and  then  printed  before  proceeding  to  the  next  time  increment.  Only 
the  quantities  required  for  the  next  calculation  are  stored. 

In  the  following  sections  of  this  report  the  procedure  is  developed  analytically,  the 
equations  used  directly  in  the  computer  program  are  indicated,  the  program  itself  is  described, 
and  the  program  is  evaluated  with  a series  of  test  problems. 

MATHEMATICAL  MODEL  FOR  SHIP  HULL 
STRUCTURAL  MODEL 

The  ship  hull  is  considered  as  a nonuniform  beam  divided  into  20  or  less  equal  sections 
or  elements.  Nodes  or  coordinates  are  used  at  the  ends  of  each  section  as  shown  in  Figure  1. 

The  mass  m.  bending  rigidity  EI/C.  and  shear  rigidity  KAG/C  of  the  sections  are  numbered  as 
shown  in  Figure  1.  The  deflections,  y^  are  taken  at  the  nodes  at  the  ends  of  the  sections. 

The  slamming  forces  F are  considered  as  discrete  forces  acting  vertically  at  the  nodes.  Both 
deflections  and  forces  are  positive  upwards.  Damping  characteristics  are  discussed  later. 

UNDAMPED  EQUATIONS  OF  MOTION 

i 

The  undamped  equations  of  motion  in  matrix  form  are: 

[Ml  { y } + i K 1 {y}  = {F}  (1) 

where  [M]  is  the  mass  matrix  (including  added  mass  due  to  water) 

[ K I is  the  stiffness  matrix 
{ y}  is  the  vector  of  vertical  accelerations 
{ y}  is  the  vector  of  vertical  deflections 
{ F}  is  the  vector  of  applied  forces 
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Structural  Model  for  Vertical  Vibration 
of  a Ship  Hull 
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The  mass  and  stiffness  matrices  are  based  on  a finite-element  beam  model  and  are  taken  from 
Reference  6.  The  mass  matrix  for  a single  finite  element  between  nodes  i and  j,  referred  to 
the  coordinate  set  {y^  0j,  y^,  0. } where  0 is  slope,  is 


156 

-228 
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- 228 
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- 382 

420 
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-138 
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228 

482 

where  m is  the  mass  of  the  element,  and  8 is  the  length  of  the  element. 


The  stiffness  matrix  for  the  element,  referred  to  the  same  coordinate  set,  is 


where 


-R 


IK] 


R 


8 + 83  \ ' 
KAG  + 12  El  / 


(3) 


(4) 


The  mass  and  stiffness  matrices  for  the  entire  beam  are  assembled  by  superposition  of  the 
element  matrices.  The  resulting  matrices  for  a beam  with  M elements  and  M + 1 nodes  will 
be  of  the  order  N = 2M  + 2. 


^ Mac  Neal.  Richard  H.,  “The  NASTRAN  Theoretical  Manual  (Level  15),"  National  Aeronautics  and  Space  Administration 
SP-221101 ).  Washington,  D.  C.  (Apr  1972). 
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DAMPING  MATRIX 


To  include  damping  in  the  equations  of  motion,  consider  a viscous  damping  matrix  |C) 
so  that  the  equations  of  motion  become 

[Ml  {y}  + [C]  {y}  + [K]  { y } = {F}  (5) 

The  normal  mode  method  requires  that  these  equations  be  uncoupled,  and  this  can  only  be 
done  if  the  damping  matrix  is  proportional  to  either  the  mass  matrix,  or  the  stiffness  matrix, 
or  a linear  combination;  see  Reference  7 or  section  about  uncoupling  equations  of  motion. 

[Cl  = a [Ml  + b [K]  (6) 

This  requirement  is  broad  enough  to  accommodate  the  types  of  damping  normally  used 
with  hull-vibration  damping.  Reference  8 considers  two  types  of  damping  coefficients  c, 

(1)  Rayleigh  damping  in  which  c/p  is  a constant  CR  where  p = m/5  is  the  mass  per  unit  length, 
and  (2)  frequency  dependent  damping  in  which  c/pw  is  a constant  CF  where  co  is  the  circular 
frequency.  Both  of  these  are  mass  proportional  and  can  be  expressed  in  terms  of  [ M ) . For 
Rayleigh  damping 


[Cl  = <Cr/6)  [Ml,  a = CR/S,  b = 0 

For  frequency  dependent  damping  it  is  more  complicated,  since  the  u varies.  In  this  case, 
we  use  the  factor  CF/ 5 but  after  uncoupling  the  modes,  each  modal  damping  constant  is 
multiplied  by  the  natural  frequency  of  that  mode;  see  section  about  uncoupling  equations  of 
motion. 

Experimental  data  (vibration  generator  tests  and  anchor  drops)  cited  in  Reference  8 
indicate  that  frequency  dependent  damping  is  more  appropriate  and  that  the  value  C,  = 0.03 
is  an  average  value  for  several  ships  tested. 

Physically  mass  proportional  damping  corresponds  to  dampers  connected  between  the 
nodes  and  an  inertial  reference. 

A portion  of  the  total  hull  damping  must  be  due  to  hysterisis  which  would  be  represented 
by  dampers  working  against  the  relative  rotational  velocities  (changes  in  slope)  and  shear  velo- 
cities of  adjacent  elements.  Although  this  type  of  damping  has  not  been  used  in  many 


^Hurty,  W.  C.  and  M.  F.  Rubinstein.  “Dynamics  of  Structures,"  Prentice-Hall,  Inc.,  Fnftlewood  Cliffs,  N.  J.  (I%4). 
SMcGoldrick,  R.  T.,  "Ship  Vibration,"  David  Taylor  Model  Basin  Report  1451  (Dec  1%0). 


calculations,  further  development  may  justify  its  use,  and  it  is  included  in  the  program.  The 
moment  and  shear  transmitted  between  two  sections  by  the  bending  rigidity  El  and  shear 
rigidity  KAG  is  dependent  on  the  slopes  and  deflections  of  adjacent  elements.  If  a damper 
also  transmits  a moment  and  shear  between  two  elements  (but  as  a function  of  the  time  rate 
of  change  of  slopes  and  deflections)  then  the  damping  matrix  must  take  the  same  form  as 
the  stiffness  matrix.  If,  in  addition,  the  values  of  the  coefficients  are  proportional  to  the 
rigidities,  then  the  damping  matrix  can  be  written 

1C]  = b IK] 

NATURAL  FREQUENCIES  AND  MODE  SHAPES 
MODES  OF  VIBRATION  OF  A FREE  FREE  BEAM 

The  first  step  in  solving  the  equations  of  motion  of  the  idealized  beam  is  the  solutioi 
of  the  free,  undamped  vibration  problem. 

[Ml  { y } + |K]  { y } = {0}  (7) 


We  seek  solutions  of  the  form 


{y} 


{ x } 


itJ  t 


(8) 


Equation  (7)  becomes 


-wn2  [M]  {X}  + IK]  {X}  = {0}  <‘> 

There  are  as  many  solutions  for  ion  as  there  are  degrees  of  freedom.  Each  value  of  u>n  is  an 
eigenvalue  (natural  frequency)  of  the  system,  and  is  dependent  only  on  the  masses  and  stiff- 
nesses of  the  system.  For  each  eigenvalue  there  is  a corresponding  vector  (X)  which, 
•ogether  with  the  eigenvalue,  satisfies  Equation  (9).  The  vector  (X)  is  an  eigenvector  and 
represents  the  mode  shape.  The  mode  shapes  can  be  determined  only  to  within  a multiplica- 
tive constant.  The  matrix  [X]  formed  by  the  column  vectors  {X}  is  called  the  eigenvector 
matrix. 

If  the  constraints  are  such  that  rigid  body  motion  is  possible,  wn  = 0 will  be  a solution 
with  a multiplicity  corresponding  to  the  number  of  rigid  body  modes.  A free-free  beam  con- 
sidered in  the  vertical  direction  has  two  rigid  body  modes,  one  in  translation  and  one  in 
rotation. 
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The  two-noded  mode  will  be  the  first  of  the  remaining  modes,  the  three-noded  mode 
will  be  the  second,  etc. 


DETERMINANT  METHOD  OF  SOLUTION 

To  solve  Equation  (9),  we  can  combine  the  mass  and  stiffness  properties  into  one 
operator  matrix. 

([K]  - wn2  [M])  {X}  = {0}  (10) 

This  represents  a set  of  N homogeneous  algebraic  equations  in  X.  For  { X } to  have  non- 
zero solutions,  the  determinant  of  the  operator  matrix  must  be  zero. 

|[K]  - u>n2  I M ] | = 0 (||) 

Trial  values  of  cun  are  substituted  into  Equation  (11),  and,  by  interpolation,  the  values  of  a>n 
that  cause  the  determinant  to  be  zero  are  found;  they  are  the  natural  circular  frequencies  of 
the  system. 

To  obtain  the  eigenvectors,  the  wn  are  substituted  back  into  Equation  (10).  and  the 
vector  { X } is  found  for  each  u>n.  If  we  let 

ID)  = <(K)  - wn2  (M|) 

then  Equation  (10)  can  be  written 

[D]  {X}  = {0} 


I 


which  when  expanded  takes  the  form 


1 1 Xl 

+ d12 
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II 
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Since  mode  shapes  can  only  be  determined  within  multiplicative  constants  we  can  specify  any 
one  Xj  in  each  mode.  We  find  it  convenient  to  let  XN  = 1 in  all  modes.  Then  the  first 
N - 1 equations  are  written 


+ d12x2 

+ . . . . 

' • • + DI,  N - 

1 XN-I 

D1N 

'21  X. 

+ d22  X2 

+ . . . , 

• • • + D2,  N - 

1 XN  - 1 

Z 

N 

Q 

1 

ii 

- 1.  1 *1  + - 1.  2 *2  + + - 1.  N - 1 XN  - 1 - 1 . N 

(13) 

These  are  solved  as  a set  of  nonhomogeneous  algebraic  equations,  and,  after  the  vectors  {X} 
are  found,  each  can  be  substituted  into  the  Ntfi  equation  to  check 

DN1  X,  + DN2  X2  + + Dn  n _ i XN  _ [ + dnn  = 0 

This  “determinant  method”  has  several  advantages  over  other  commonly  used  methods 
for  finding  eigenvalues  and  eigenvectors.  Possibly  the  greatest  advantage  is  that  it  is  not 
necessary  to  reduce  the  order  of  the  dynamic  matrix  by  the  number  of  rigid  body  modes, 
thereby  obtaining  the  equations  of  motion  in  generalized  coordinates.  While  this  may  not  be 
difficult  for  the  beam,  if  the  program  is  later  adapted  to  other  transient  problems,  each  type 
of  constraint  must  be  treated  separately. 

Other  advantages  are  that  it  is  more  accurate  for  high  modes  than  some  other  methods 
and  that  only  the  modes  desired  need  be  calculated. 

CALCULATION  OF  RESPONSE 


ORTHOGONALITY 

In  this  section  we  will  show  that  the  eigenvectors  are  orthogonal  with  respect  to  the 
matrices  [Ml  and  [K] . 

If  we  write  Equation  (9)  for  the  ith  mode,  we  have 

cu2  [Ml  {X(i>|  = [K]  {X(i,|  (14) 
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where  {X(1)}  is  the  mode  shape  (eigenvector)  for  the  ith  mode.  Premultiplying  this  by 

7 

{ X(J) } , the  transposed  eigenvector  for  a different  mode,  we  get 

w2  {X«>}T  [M]  {X(i)}  = (X<J>}T  [K]  {X(i)}  (15) 

Next,  take  the  transpose  of  both  sides  of  Equation  (15)  and  apply  the  rule  which  states  that 
the  transpose  of  a matrix  product  is  equal  to  the  product  of  the  transposed  matrices  in 
reverse  order. 


<o2  (XV  IM)  IX<»}  = {X«>}T  IK)  {X0)}  (16) 

We  have  used  the  fact  that  (M)  and  )K)  are  symmetric  and  therefore  equal  to  their  trans- 
poses. Now  we  write  an  equation  identical  to  Equation  (15)  except  interchanging  indices. 

co,2  {Xt0}T  [M J {X<J>}  = {X(i,}T  (K)  {X01}  (17) 

Subtract  Equation  (17)  from  Equation  (16) 

(co2  -w/)  {X(i)}T  [Ml  {X<>>}  = {0}  (18) 

For  cOj  =£ 


(X(i)}T  )M)  { X<j) } = { 0 } (19) 

Also,  from  Equation  (17)  we  have 

{X°'}T  |K|  {X^1}  = {0}  (20) 

Equations  (19)  and  (20)  reflect  the  definition  of  orthogonality  with  respect  to  a weighting 
matrix,  in  this  case  ( M ] or  (K).  Note  that  if  two  natural  frequencies  are  the  same,  their 
mode  shapes  are  not  necessarily  orthogonal. 


NORMALIZATION 

We  mentioned  before  that  an  eigenvector  represents  a characteristic  pattern  of  relative 
displacements  associated  with  a particular  mode  of  vibration.  This  means  that  we  can  multi- 
ply all  the  components  of  an  eigenvector  by  any  constant,  and  it  will  still  represent  the  same 
mode  shape.  We  are  able  to  use  this  fact  to  simplify  our  analysis. 

T 

Premultiplying  Equation  (14)  by  {X(l)}  , the  transposed  eigenvector  for  the  ith  mode, 
we  obtain 


u;2  {X(i)}T  [M]  {X(i)}  = {X(i)}T  [K]  { X°>}  (21) 

r 

First  we  will  consider  the  triple  matrix  product  on  the  left  side.  It  can  be  shown  that  the 
product  will  be  some  constant,  say  ms 


{X(i)}T  [M]  {X"'}  = m,  (22) 

We  wish  to  normalize  the  vector  { X(i) } and  its  transpose  by  multiplying  it  by  a constant  n; 
so  that  ntj  becomes  unity 


n2  {X(0}T  [M]  {X(,)}  = 1 


(23) 


Written  in  subscript  notation 


ni2  Z Z "V  Xii  Xk.  = 


(24) 


i k 


where  m.  is  an  element  of  [Ml.  and  XH  is  the  j ih  component  of  {X(i)}.  Solving  for  n, 

JK  J*  • 


n- 


|/XT 


(25) 


m. 


j k 


jk 


Xji  xki 


Let  the  normalized  eigenvector  be  denoted  by  { X,,'1  } , so  that 


{ Xn(i) } = n(  { X<‘> } 


(26) 


If  we  write  Equation  (21)  using  our  normalized  eigenvectors,  we  get 


W2  {V>}T  lMl  {Xn<0>  = tXn(°}T  IK)  tXn(°>  <27> 


We  have  shown  that 


{X"}T  [M]  {Xn(i)}  = 1 


(28) 


which  means  * 


(V1'}1  IK1  Un(i)}  = w2  (29) 

If  we  assemble  the  normalized  eigenvector  matrix  [X^l,  we  can  show  by  use  of  Equations 
(19),  (20).  (28),  and  (29)  that 


[XJT  [M]  [XJ  = [I]  (30) 

lXn]T  IK1  [Xn]  = (co2)  (31) 

where  |u>2)  is  a diagonal  matrix  of  the  squares  of  The  last  two  equations  are  useful  in 
our  normal  mode  analysis. 

The  eigenvector  matrix  for  the  free-free  beam  is  obtained  by  assembling  the  (N  - 2) 
and  the  rigid  body  mode  shapes  into  an  (N  x N)  matrix. 

UNCOUPLING  THE  EQUATIONS  OF  MOTION 

Our  next  task  is  to  uncouple  the  equations  of  motion  so  that  we  have  a series  of  equa- 
tions (one  for  each  mode)  resembling  that  of  a simple  oscillator.  We  start  with  Equation  (5) 

[M]  {y } + [Cl  {y  } + [K]  (y  } = {F}  (5) 

where  {y  } and  its  derivatives  and  { F } are  functions  of  time.  The  displacement  response 
may  be  expressed  as  a superposition  of  normal  modes 

{ y(t) } = [Xn]  { q(t) } (32) 
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where  {q}  is  a vector  of  time-dependent  generalized  coordinates.  They  are  called  normal 
coordinates  since  they  are  associated  with  the  normalized  modal  matrix  [XJ.  Substitution 
gives 


IM)  [XJ  {q}  + 1C]  [XJ  { q } + IK)  [XJ  (q } = {F}  (33) 

Premultiply  this  equation  by  IXJT 

IX„)T  IM)  [XJ  {q}+  |XJT  [Cl  [XJ  {q  }+  { X„  } T [KJ  [XJ  {q  } = [XJT  (F> 

(34) 


Applying  Equations  (30),  (31),  and  (6)  yields 

[I)  {q  } + a [I  J { q } + b [or ) {q}+[ur)  {q}  = [Xn)T  {F>  (35) 


If  we  let 


(P(t)}  = [XnlT  {F(t»  (3b) 

1 G ) = a [1]  + b |w:]  (37) 

Equation  (35)  becomes 

(q}+  IG)  {q}  + [co2]  {q}  = { P( t ) } (38) 

The  matrices  [ G J and  [co2)  are  diagonal  so  that  Equation  (38)  can  be  written  as  a system  of 
uncoupled  equations 


qj(t)  + Gjq/t)  + cu,2  qj(t)  = P,(t),  i = 1,  2 N (39) 

where 

Gj  = a + bu>2  (40) 

Equations  (39)  have  the  form  of  the  equations  of  motion  of  damped  simple  oscillators, 
and  each  equation  represents  a particular  mode  of  vibration.  The  solution  of  Equations  (39) 
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is  well  known,  and  there  are  several  methods  that  could  be  used.  Our  computer  program 
uses  a time-marching  method. 

TIME-MARCHING  METHOD 

There  are  several  time-maching  methods  which  approximate  the  differential  equations  of 
motion  with  finite-difference  equations  so  that  the  displacement  at  a particular  time  t , is 
an  algebraic  function  of  the  displacement  at  time  tn,  the  parameters  of  the  system,  and  the 
applied  force.  These  methods  appear  simpler  than  transform  methods  which  involve  a large 
amount  of  numerical  integration. 

The  particular  time-marching  method  used  herein  is  one  developed  by  Chan,  Cox,  and 
Benfield.9  In  this  method  we  consider  time  to  be  divided  into  many  small,  equal  time 
increments,  h,  where  h = tn  + 1 - t„  = tn  - tn  ,.  The  modal  displacements  qn+,  at  the  time 
tn  + 1 are  calculated,  then  transformed  back  to  physical  coordinates. 

The  equations  are  derived  in  Reference  9;  only  the  results  are  stated  here.  The  main 
recurrence  formula  for  the  solution  of  Equation  (39)  is 

Zqn+1  = Aqn  - Wq„  , + 0h2P„  + 1 + <1  - 2 0)  h2Pn  + 0h2Pn  , , n > 0 (41) 


where 


Z = I + (h/2)C,  + 0h2  w2 

A = 2 (1  20)  h2  co2  (42) 

W = 1 (h/2)G  + 0h2  oj2 

Equation  (41)  must  be  applied  to  each  of  the  modal  displacements;  however,  modal  subscripts 
have  been  omitted  to  avoid  confusion  with  the  time  subscripts.  The  quantity  0 can  be  assigned 
any  value  between  0 and  1/4.  Different  values  of  0 correspond  to  different  approximating 
assumptions  about  the  acceleration  between  time  steps  (such  as  a constant  acceleration  within 
each  increment  or  a linear  function).  The  effect  of  0 on  accuracy  is  discussed  in  the 
Evaluation. 

Equations  (41 ) and  (42)  are  not  valid  for  the  first  time  increment.  The  following 
equation,  which  includes  initial  conditions  must  be  used. 


9Chan,  S.  cl  al„  "Transient  Analy  sis  of  Forced  Vibrations  of  Complex  Structural-Mechanical  Systems.”  Journal  of  the 
Royal  Aeronautical  Society,  Vol.  66.  No.  457  (Jul  1962). 
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r 


Zq*  = Qq0  + Rq„  + SP0  + 0h2P, 


where 


Q = I + hG/2  - (1/2  - 0)  h2  w2  - (1/4  - 0)  h3  Gw2 
R = h - (1/4  - 0)  h3  G2 
S = (1/2  - 0)  h2  + (1/4  - 0)  h3G 


(43) 


(44) 


q0  = initial  modal  displacement 
q0  = initial  modal  velocity 

After  modal  displacements  at  a particular  time  have  been  determined  by  the  time- 
marching equations,  the  nodal  displacements  are  obtained  for  Equation  (32). 

y =[Xnl  q (32) 


From  the  nodal  displacements,  we  can  obtain  the  velocities  and  accelerations.  The  velocity 
is  found  as  a second  order  approximation  of  the  slope  of  the  displacement,  and  the  accelera- 
tion is  calculated  from  the  slope  of  the  velocity  or  curvature  of  displacement 


n + 1 


*^,x" 


4X_  + 


3 V,  > 


(45) 


n + 1 


(Xn_,  -2Xn  + Xn+1) 


The  moment  in  the  element  between  nodes  i and  j is  approximated  by 


El. 


m,  * T <9,  - *.» 


The  stress  Sb  due  to  bending  in  the  element  between  nodes  i and  j is  given  by 

% 


Vj  Mij  c^ij 


(4M 


(47) 


(48) 
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where  c is  the  vertical  distance  from  the  neutral  axis  to  the  point  where  the  stress  is  being 
calculated. 

The  shear  force  V at  the  ]th  node  is 


M., 


(49) 


The  shear  stress  at  any  station  near  the  neutral  bending  axis  is 


(50) 


where  KA  is  the  effective  area  at  the  station.  (In  our  model  the  effective  areas  are  given  by 
elements,  so  the  KA  of  either  adjacent  element  could  be  used.) 


DESCRIPTION  OF  PROGRAM 
SUMMARY  OF  APPLICABLE  EQUATIONS 

Before  discussing  the  details  of  the  program,  this  section  lists  the  equations  used  in  the 
program  trom  the  preceding  analysis  and  indicates  briefly  what  steps  are  necessary. 

first  the  mass  and  stiftness  matrices  are  obtained  from  the  superposition  of  the  individual 
element  matrices.  Equations  (2)  and  (3). 


156 

— 22C 

54 

13C 

m£ 

— 22C 

4C2 

- 13C 

- 3S2 

420 

54 

- 1 3P 

156 

1 36 

- 3C2 

22C 

4C2 
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Then  solutions  lire  obtained  lor  natural  frequencies  con  using  Equation  (II) 

UK)  wn2|M)|=0  (M) 

Tlie  eigenvectors  can  then  be  found  using  liquation  (10) 

( I K | con:  | M | ) { X } = 0 ,10) 

where  XN  = I . 

Then  the  eigenvector  matrix  must  be  normalized  by  means  of  the  following  two 
eq  uations 


{ Xn“> } = n.  {X"'} 
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The  nodal  exciting  forces  must  be  transformed  into  modal  exciting  forces 


{P>  = [X„]T  {F} 


(36) 


Tlie  modal  damping  is  calculated  from 


G,  = a + bur  (40) 

At  this  point  we  have  all  we  need  to  start  our  time-marching  procedure.  First  we  note 
that  the  starting  Equation  (43)  and  the  main  recurrence  Equation  (41)  require  the  calculation 
of  the  following  parameters 


Z * I + (h/2)G  + Ph2  or 

A = 2 - (1  - 20)  h2  ur  (42) 

W = 1 - <h/2)G  + 0h2  co2 


Q = 1 + hG/2  - (1/2  - (3)  h2  w2  - (1/4  - p)  h3  w2G 
R = h - (1/4  - p)  h3  G2  ,44) 

S = (1/2  - p)  h2  + (1/4  - p)  h3G 

Although  we  have  omitted  modal  subscripts,  these  parameters  must  be  calculated  for  each 
mode. 

Next  we  start  a series  of  calculations  which  must  be  performed  at  each  time  increment. 
The  first  step  is  either  the  starting  equation 

Zqj  = Qq0  + Rq0  + SPo  + &'1?\  (43) 

or  the  main  recurrence  formula 

= - wqn  t +0h2Pnt|  +(l  2P)  h2P„  + h2Pn_,  (41) 
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Alter  tlio  moil  it  I displacements  at  time  tnM  have  been  calculated,  we  calculate  the  physical 
or  nodal  displacements  lor  time  tnM  using 

{y}=|XJ{q}  ,32, 


rhe  nodal  velocity  and  acceleration  are  calculated  from 


X 


n + l 


4*n  + .<xnt|) 


(45) 


Xntl  = ( l/lr ) ( X 


n l 


2X„ 


+ x,ltl> 


(4(i) 


T he  bending  moments  are  calculated  from 


= T 


flic  bending  stresses  are  calculated  from 


(47) 


(4K) 


The  'hear  force  is  cal<  ulated  from  liquation  <4‘)>  and  the  shear  stress  from  liquation  (50) 
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j 


M , M 


ik 
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( 4° ) 
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(50) 


After  these  calculations  the  output  for  each  desired  station  for  each  time  increment  is 
printed,  and  the  program  begins  the  loop  again  for  the  next  time  increment. 


MAIN  PROGRAM 

Appendix  A is  a flow  chart  of  the  program  SLAM  and  indicates  how  the  equations  of 
the  previous  section  are  incorporated  into  the  program.  SLAM  is  designed  to  calculate  the 
vertical  transient  response  of  a beam-type  structure  with  up  to  20  elements.  Vertical  forces 


can  be  applied  to  the  first  10  stations.  At  each  of  the  stations  the  force  can  be  specified  by 
as  many  as  100  points  with  respect  to  time.  The  response  is  calculated  for  every  hundredth 
of  a second  which  is  adequate  to  describe  motions  with  frequencies  at  least  up  to  10  hertz. 
Output  for  every  increment,  or  for  every  NPRINT  increments,  can  be  printed  as  desired.  The 
user  specifies  the  number  of  modes  to  be  included  in  the  response  calculation  and  can  opt  to 
have  the  natural  frequencies  and  mode  shapes  printed  out.  At  each  station  for  which  output 
is  requested,  the  displacement,  acceleration,  bending  moment,  shear  force,  and  bending  and 
shear  stresses  are  printed  out.  Since  the  output  is  printed  by  stations  and  the  bending  moment, 
bending  stress,  and  shear  stress  are  associated  with  the  elements,  those  quantities  are  given  for 
the  element  just  forward  and  just  aft  of  the  station  named.  SLAM  will  account  for  any  of 
three  types  of  damping: 

1.  Mass  proportional 

2.  Stiffness  proportional 

3.  Mass  and  frequency  proportional 

For  the  input,  any  consistent  set  of  units  can  be  used,  except  that  the  value  of  E used  can 
be  a different  set  of  units.  The  stresses  calculated  will  be  in  that  set  of  units. 

SUBROUTINES 

FORCES  - This  subroutine  takes  the  forces  as  plotted  for  the  program  input,  inter- 
polates where  necessary,  and  calculates  the  modal  forces  which  are  used  in  the  main  program. 

MATINS  This  subroutine  is  in  the  Center  library  of  subroutines  and  can  either  find  the 
determinant  of  a square  matrix  |A|  or  indicate  if  it  is  singular  and  either  find  the  inverse  of 
the  matrix  [A]  or  solve  a set  of  simultaneous  linear  equations  of  the  form  |A]  {X}  = |B). 

In  SLAM  it  is  used  first  to  find  determinants  in  the  natural  frequency  part  of  the  program. 

Then  it  is  used  to  solve  simultaneous  equations  in  obtaining  the  mode  shapes. 

INPUT/OUTPUT 

Sample  input  for  SLAM  is  given  in  Figure  2.  The  necessary  format  and  description  of 
the  input  is  given  in  Table  2. 

Sample  output  for  SLAM  is  given  in  Figure  3. 

EVALUATION 

ANALYTICAL  APPROACH 

The  analytics  used  in  SLAM  were  selected  specifically  for  the  ship-slamming  applica- 
tion; therefore,  they  should  provide  maximum  flexibility  and  utility. 
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Figure  2 - Sample  Input  for  SLAM  (See  Table  2) 
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TABLfc  2 KXPLANATION  OK  SLAM  INPUT 


Line 

Numbei 

Format 

Description 

1 

(12.20(1  X . 1 2 ) 

N - Number  of  sections  (20)* 

NMODE  Number  of  modes  to  be  calculated  (101* 

lOPT  (11  - For  input 

•OPT (2)  For  natural  frequencies 

IOPTI3)  For  mode  shapes 

IOPT (4 ) For  transient  response 

( 1 OP T ( 1 1 indicates  what  output  is  desired-1  if  desired.  0 if 
not  1 

NPTS  Number  of  points  .n  the  force  versus  time  plot  including 

t 0(1001* 

NFOR  Number  of  sections  to  which  force  is  applied  (101* 

NOUT  Numbei  of  sections  for  which  output  is  desired  (201* 

NPRINT  Number  of  calculations  for  each  printout 

2 

(8E10  3) 

DELX  Length  of  hull  sections 

DAMPM  Mass  proportional  damping  constant 

DAMPK  Stiffness  proportional  damping  constant 

DAMPF  Mass  and  frequency  proportional  damping  constant 

TSTOP  Time  of  last  desired  response 

E Modulus  of  elasticity 

S Shear  modulus 

3 22 

(86 10  3) 

EMASSII)  Mass  of  hull  elements 

EE llll  El  of  hull  elements 

EKACi'fi  KAG  of  hull  elements 

D! ST ( 1 1 Distance  from  neutral  bending  axis  to  desired  stress 

location 

(Consists  of  N lines  starting  at  bow  ) 

23  93 

(8E10  3) 

TPT(I)  - Time  of  each  force  versus  time  promt 

FPT(I.J)  Force  on  the  NFOR  sections 

(Consists  o*  NPTS  lines,  each  one  gives  the  time  and 

NFOR  forces  applied  at  that  time.  Includes  t - 0 at 
which  time  all  forces  must  equal  zero,  also,  must  give 
forces  up  to  t TSTOP. 

94 

(I2.20MX.I2) 

LOUTII)  Station  numbers  for  which  output  is  desired 

(Consists  of  NOUT  station  numbers.  Do  not  include  0 
and  N which  printout  automatically.) 

•Indicates  maximum  number  allowed 

Figure  3 - Sample  Output  for  SLAM 
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The  mathematical  model  selected  is  a finite-element  representation  of  a beam  with  both 
bending  and  shear  rigidities,  which  should  give  very  accurate  results.  A matrix  approach  is 
used  so  that  more  complicated  hull  forms  (such  as  catamarans  or  surface-effects  ships)  will 
require  only  reprograming  the  matrix  assembly  portion  to  relied  the  more  complex  structural 
model.  Once  the  mass  and  stiffness  matrices  are  formed  most  of  the  rest  of  the  program  can 
be  used  unchanged. 

The  determinant  method  for  eigenvalue  extraction  has  been  selected  because  it  is  accu- 
rate for  higher  as  well  as  lower  modes  and  because  with  this  method  the  rigid  body  modes 
of  the  beam  do  not  necessitate  preliminary  manipulation  as  some  methods  do. 

The  normal  mode  method  allows  the  user  to  calculate  only  the  modes  of  vibration  of 
interest.  Rigid  body  motion  is  automatically  eliminated  If  desired,  the  response  in  each 
mode  can  be  printed  out  separately  with  a minor  modification. 

The  time-marching  method  makes  it  possible  to  conveniently  input  a detailed  description 
of  the  slamming  forces  with  respect  to  time  and  length  along  the  ship.  The  parameter  0 can 
be  varied  to  minimize  errors  in  the  maximum  amplitude  or  in  the  period.  Table  3 (developed 
from  Reference  10)  indicates  the  tradeoff  that  has  to  be  made.  SI  AM  uses  0 = I 8.  Table  4 
(also  from  Reference  10)  indicates  the  highest  frequency  for  a stable  and  converging  solution. 

In  the  ship  slamming  problem  the  participation  of  the  higher  modes  has  many  uncer- 
tainties associated  with  it.  first  the  calculated  slamming  forces  are  normally  derived  from  a 
statistical  approach.  This  means  that,  although  a pressure  is  generally  triangular  with  respect 
to  time  and  is  represented  as  such,  a particular  slam  may  be  closer  to  a half  sine,  or  may  be 
skewed  to  the  left  or  right,  or  may  have  some  other  shape.  This  changes  considerably  the 
excitation  at  the  higher  frequencies,  both  in  magnitude  and  in  phase.  The  maximum  response 
of  the  ship  will  depend  on  how  the  higher  modes  respond  relative  to  the  larger  fundamental 
response,  l or  example,  if  both  modal  responses  peak  together,  the  maximum  response  would 
be  significantly  higher  than  if  the  modal  responses  subtracted  from  each  other.  Another 
complication  is  that,  at  least  in  vibration  generator  tests  on  higher  modes  of  ships,  the  meas- 
ured mode  shapes  are  sometimes  significantly  different  from  the  calculated  normal  mode 
shapes.  Additional  measured  data,  compared  with  calculations,  should  help  to  describe  the 
effect  of  the  higher  modes. 

TEST  PROBLEMS 

To  test  the  accuracy  ot  SI  AM  and  to  demonstrate  its  use  and  characteristics,  a series  of 
test  problems  was  solved. 


^Ncwmark.  V M.  "Compulation  <»t  iKnamu  Structural  Response  in  the  Ranee  Approaching  I ailmo.  ' Symposium  on 
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TABLE  3 EFFECT  OF  FREQUENCY  ON  ERRORS  DUE 
TO  NUMERICAL  PROCEDURE  IN  SLAM 


Hz 

o 

II 

*3  1/12 

(3-V8 

<3  - 1/6 

<3  = 1/4 

RELATIVE  ERRORS  IN  MAXIMUM  RESPONSE 
TO  AN  INITIAL  VELOCITY 


BB 

0 008 

1 

0 034 

EH 

0.166 

0.614 

0 306 

0 004 

0 

0017 

0 

0 073 

0 

0 122 

0 

RELATIVE  ERRORS  IN  PERIOD 


-0.004 

-0.017 


-0  0001 

0.002 

0 004 

-0  0003 

0 008 

0.017 

-0.006 

0 028 

0.059 

-0015 

0.038 

0087 

TABLE  4 STABILITY  AND  CONVERGENCE  LIMITS  FOR 
P METHOD  IN  SLAM 


Stability  Limit  in  hertz 


Convergence  Limit  in  hertz 


*3  = 0 *3  = 1/12 


31  8 


*3  = 1/4 


45.0 

55.1 

oo 

45  0 

38.9 

31.8 

The  first  series  of  problems  involved  a 600-foot  uniform  steel  beam  weighing  1000  tons 
with  a bending  rigidity,  El,  of  1010  ton-ft2.  The  shear  rigidity  was  assumed  to  be  I021 
(very  larger.  The  end  of  the  beam  was  excited  laterally  by  a 10,000-ton  step  function,  and 
the  undamped  transient  response  for  the  first  two  modes  was  calculated.  The  exact  solution 
for  displacement  and  acceleration  of  the  end  of  the  continuous  beam,  obtained  from  Refer- 
ence 7,  is  given  by  the  solid  lines  of  Figure  4.  The  beam  was  then  divided  into  20  sections, 
and  solutions  were  obtained  from  SLAM  which  are  given  by  the  circles  in  Figure  4. 

In  addition  to  the  uniform  beam,  a MARINI  R-('lass  hull  was  used  for  some  test  prob- 
lems. Figures  5 through  8 show  the  response  of  the  MARINER  hull  to  slamming  forces 
calculated  using  SLAM  for  a slightly  more  severe  than  State  6 sea.  Frequency  and  mass 
dependent  damping  was  used,  causing  the  higher  mode  components  to  decay  before  the 
fundamental  component  as  shown  from  experimental  results.  Figures  5 and  6 show  the 
response,  considering  only  the  first  three  modes.  Figures  7 and  8 shew  the  response,  con- 
sidering the  first  10  modes.  The  displacement  and  accelerations  were  plotted  for  the  bow  and 
Station  7.  The  bending  moment  and  shear  forces  were  plotted  at  Station  7 and  amidships. 

CONCLUSIONS 

Normally  a ship  hull  will  be  divided  into  20  sections  for  slamming  calculations  for  two 
reasons:  the  parameters  are  likely  to  be  listed  at  20  stations  for  other  design  studies,  and  20 
stations  is  an  appropriate  number  for  describing  the  forces  applied  to  the  hull.  The  accuracy 
of  the  solution  of  the  uniform  beam  (Figure  4)  indicates  that  numerical  procedures,  using  a 
20-element  representation,  are  sufficiently  accurate  for  any  projected  use  in  connection  with 
conventional  ship  slamming  involving  the  lower  modes.  A 10-element  model  would  probably 
do  almost  as  well  in  the  first  three  modes;  however,  the  force  breakdown  may  be  too  coarse. 

The  user  should  study  Figures  5 through  8 to  determine  the  effect  of  including  higher 
modes.  The  displacement  and  bending  moment  do  not  reflect  the  higher  modes  as  much  as 
the  acceleration  and  shear  force.  The  number  of  modes  requested  should  be  determined  by 
which  parameters  are  important  for  that  particular  application  and  by  the  accuracy  needed. 
Not  much  is  gained,  however,  by  requesting  more  modes  than  half  the  number  of  elements. 

The  flexibility  of  the  analytics  used  in  SLAM  allows  the  user  to  study  intermediate 
steps  in  the  structural  solution:  natural  frequencies,  mode  shapes,  and  individual  modal 
responses.  SLAM  is  flexible  also  in  that  it  can  be  easily  adapted  to  handle  structural  con- 
figurations other  than  a beam. 
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Bending  Moment  and  Shear  Force  of  a MARINER-Class  Hull.  Considering  Three  Modes 
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Figure  7 Displacement  and  Acceleration  Response  of  a MARINER-Class  Hull.  Considering  10  Modes 


APPENDIX  A 


FLOW  CHART  OF  SLAM 
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